Analyse des données RNAseq DepMap (CCLE ; cancer) et des 57 epigenomes de Roadmap
# W O R K I N G D I R E C T O R Y
main.dir = "/shared/projects/chrom_enhancer_bc"
setwd(main.dir)
data.dir = file.path(main.dir,"data/expression/rdata")
# I M P O R T
library(pheatmap )
library(ggpubr)
library(sva)
library(edgeR)
library(DESeq2)
library(DT)
library(stringr)
library(dplyr)
library(reshape2)
library(FactoMineR)
library(factoextra)
library(grid)
source("./scripts/Normalization.R")
load(file.path(data.dir,"cancer_healthy_expression.RData"))
# F U N C T I O N S
# Centering by the median
MedianCentering <- function(x){
(x - median(x))
}
get_color_scale_pheatmap = function(matrix,ncolors = 50,threshold = 1){
range <- max(abs(matrix))
myBreaks = seq(-range, range, length.out = ncolors)
myBreaks[myBreaks > -threshold & myBreaks < threshold ] = 0
myBreaks = unique(myBreaks)
paletteLength <- length(myBreaks)
myColor <- colorRampPalette(c("blue", "white", "red"))(paletteLength-1)
return(list("colorscale" = myColor, "breaks" = myBreaks)) }
# D A T A P R O C E S S
# change format
cell_lines = names(expr)[! names(expr) %in% c("gene_name","ENSEMBL_id")]
expr = select(expr, -gene_name)
# as numeric
expr <- cbind("ENSEMBL_id"=expr$ENSEMBL_id,
mutate_all(select(expr,all_of(cell_lines)), function(x) as.numeric(as.character(x))))
# Remove lowly expressed genes + 1 duplicated id
expr <- expr %>%
filter(! duplicated(ENSEMBL_id))
row.names(expr) = expr$ENSEMBL_id
expr$ENSEMBL_id = NULL
Row_sums = rowSums(expr)
filtered_expr = expr[Row_sums > 3*ncol(expr),]
DT::datatable(head(filtered_expr))
# Parameters
design = c(rep("CCLE", 20 ), rep("RoadMap",length(cell_lines)-20))
count.matrix = data.matrix(filtered_expr)
All_normalisations = function(count.matrix,design){
norm_list = list()
tools=c("TMM", "TMMwsp", "RLE", "upperquartile")
# edgeR
for (tool in tools){
message(tool)
edgeR.dgelist = DGEList(counts = count.matrix, group = factor(design))
nf = calcNormFactors(edgeR.dgelist, method=tool)
edgeR_norm <- cpm(nf, normalized.lib.sizes=TRUE, log = TRUE)
norm_list[[tool]] = edgeR_norm
}
#DESeq2
message("DESeq2")
deseq_norm = tools.norm.RNAseq(count.matrix, tool = "vst2", design = design )
norm_list[["DESeq2"]] = deseq_norm
return(norm_list)
}
norm_list = All_normalisations(count.matrix,design)
## TMM
## TMMwsp
## RLE
## upperquartile
## DESeq2
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Applying limma batch correction
norm_list_limma = lapply(norm_list, function(x) removeBatchEffect(x, design) )
norm_list_limma_melted = lapply(norm_list_limma, function(x) reshape2::melt(x) )
# ComBat batch corrections : needs raw counts
load(file.path(data.dir,"ComBat_seq_correction.RData"))
norm_list_ComBat = All_normalisations(adjusted,design)
## TMM
## TMMwsp
## RLE
## upperquartile
## DESeq2
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
norm_list_ComBat_melted = lapply(norm_list_ComBat, function(x) reshape2::melt(x) )
Scatterplot Limma batch correction vs ComBat_seq
k562 = colnames(norm_list_limma[['TMM']])[grepl("K562",colnames(norm_list_limma[['TMM']]))]
Limma_k562 = lapply(norm_list_limma, function(x) select(data.frame(x), all_of(k562) ))
ComBat_k562 = lapply(norm_list_ComBat, function(x) select(data.frame(x), all_of(k562) ))
ggplot(select(data.frame(log2(count.matrix+1)), K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE,K562), aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 20, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 17.5, aes(label = ..rr.label..)) +
ggtitle("Raw counts (log2+1)")
## `geom_smooth()` using formula 'y ~ x'
#K562
ggplot(Limma_k562$TMM, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("TMM normalisation - Limma batch correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(Limma_k562$TMMwsp, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("TMMwsp normalisation - Limma batch correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(Limma_k562$RLE, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("RLE normalisation - Limma batch correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(Limma_k562$upperquartile, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("Upperquartile normalisation - Limma batch correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(Limma_k562$DESeq2, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 20, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 15, aes(label = ..rr.label..)) +
ggtitle("vst normalisation (DESeq2) - Limma batch correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(ComBat_k562$TMM, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("TMM normalisation - Combat_seq correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(ComBat_k562$TMMwsp, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("TMMwsp normalisation - Combat_seq correction ")
## `geom_smooth()` using formula 'y ~ x'
ggplot(ComBat_k562$RLE, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("RLE normalisation - Combat_seq correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(ComBat_k562$upperquartile, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 15, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 10, aes(label = ..rr.label..)) +
ggtitle("Upperquartile normalisation - Combat_seq correction")
## `geom_smooth()` using formula 'y ~ x'
ggplot(ComBat_k562$DESeq2, aes(x=K562_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE, y=K562)) +
geom_point() +
geom_smooth(method = "lm") +
stat_regline_equation(label.y = 20, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 18, aes(label = ..rr.label..)) +
ggtitle("vst normalisation (DESeq2) - Combat_seq correction")
## `geom_smooth()` using formula 'y ~ x'
Test des normalisations seules vs log2(raw+1)
ctcfl = "ENSG00000124092"
ctcf = "ENSG00000102974"
pou5f1 ="ENSG00000204531"
myc = "ENSG00000136997"
polr3gl="ENSG00000121851"
polr3g ="ENSG00000113356"
CCLE = cell_lines[1:20]
ROADMAP = cell_lines[21:length(names(expr))]
#melted_expr = reshape2::melt(edgeR_norm)
condition = ifelse(norm_list_ComBat_melted$TMM$Var2 %in% CCLE, yes = "CCLE", no = "RoadMap")
# cancer lines of interest :
of_interest = c("MDAMB231","MDAMB436","HCC1937", "HCC1806",
"HCC1954","MDAMB468","K562")
color_cell = ifelse(grepl(paste(of_interest,collapse = "|"), cell_lines ), yes = "red", no ="black")
expression_plot = function(melted_df){
plot = ggplot(data = melted_df, aes(x=Var2, y=value, color = condition)) +
geom_boxplot() +
geom_text(aes(x=Var2,
label=ifelse(grepl(polr3g,Var1),".",''))
, size= 10, color = "black") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = color_cell),
axis.text=element_text(size=6)) +
xlab("cell_lines") + ylab("normalized_expression")
return(plot)
}
compute_plot = function(liste, title = "Normalisation"){
first_part_title = "Expression boxplot :"
c = 1
plot_list = list()
for (df in liste) {
tmp_title = paste(first_part_title, names(liste)[c] ,title)
plot = expression_plot(df) + ggtitle(tmp_title)
plot_list[[ names(liste)[c] ]] = plot
c = c+1
}
return(plot_list) }
# Raw + log2
melt_raw = reshape2::melt(log2(count.matrix+1))
raw_log = expression_plot(melt_raw) + ggtitle("Expression boxplot : Raw expression (log2+1)")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
norm_list_melted = lapply(norm_list, function(x) reshape2::melt(x) )
norm_plot = compute_plot(norm_list_melted)
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
print(norm_plot$DESeq2)
raw_log
for (plot in norm_plot){
print(plot)
}
Test des batchs correction seules vs DESeq2 norm vs log2(raw+1)
# RAW
raw_log
# Norm + batch(limma)
batch_norm_melt = norm_list_limma_melted$DESeq2
expression_plot(batch_norm_melt) + ggtitle("Expression boxplot : Normalized expression + batch correction (limma)")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
# Norm + batch(ComBat)
batch_norm_combat_melt = norm_list_ComBat_melted$DESeq2
expression_plot(batch_norm_combat_melt) + ggtitle("Expression boxplot : Normalized expression + batch correction (Combat)")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
# Raw + batch(ComBat)
batch_raw_combat_melt = reshape2::melt(log2(adjusted+1))
expression_plot(batch_raw_combat_melt) + ggtitle("Expression boxplot : Raw expression (log2+1) + batch correction (Combat)")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
# raw + batch(limma)
batch_raw = removeBatchEffect(log2(count.matrix+1), design)
batch_raw_melt = reshape2::melt(batch_raw)
expression_plot(batch_raw_melt) + ggtitle("Expression boxplot : Raw expression (log2+1) + batch correction (limma)")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
Vérification avec PCA
Raw
raw_matrix = data.frame(log2(count.matrix+1))
IQRs = data.frame("IQR" = apply(raw_matrix, 1, IQR))
IQRs$IQR = as.numeric(IQRs$IQR)
pca_subset = raw_matrix[rownames(filter(IQRs, IQRs$IQR>7)) , ]
hist(IQRs$IQR)
of_interest = c("MDAMB231","MDAMB436","HCC1937", "HCC1806",
"HCC1954","MDAMB468")
design_PCA = ifelse(grepl("K562", cell_lines), yes = "K562", no = design)
design_PCA = ifelse(grepl(paste(of_interest,collapse = "|"), cell_lines), yes = "TNBC", no = design_PCA)
res.pca = PCA(t(pca_subset), graph = F)
fviz_pca_ind(res.pca,
col.ind = design_PCA,
repel = TRUE,
palette = "Set1",
ggrepel.max.overlaps = 1) + labs(title ="PCA raw matrix",)
Normalisation
norm_pca = data.frame(t(norm_list$DESeq2[rownames(filter(IQRs, IQRs$IQR>7)) , ]))
res.pca = PCA(norm_pca, graph = F)
fviz_pca_ind(res.pca,
col.ind = design_PCA,
repel = TRUE,
palette = "Set1",
ggrepel.max.overlaps = 1) + labs(title ="PCA normalisad count (DESeq2)",)
Normalisation + Batch correction (ComBat_seq)
norm_pca_combat = data.frame(t(norm_list_ComBat$DESeq2[rownames(filter(IQRs, IQRs$IQR>7)) , ]))
res.pca = PCA(norm_pca_combat, graph = F)
fviz_pca_ind(res.pca,
col.ind = design_PCA,
repel = TRUE,
palette = "Set1",
ggrepel.max.overlaps = 1) + labs(title ="PCA normalisad count (DESeq2)",)
Normalisation + Batch correction (limma)
norm_pca_limma = data.frame(t(norm_list_limma$DESeq2[rownames(filter(IQRs, IQRs$IQR>7)) , ]))
res.pca = PCA(norm_pca_limma, graph = F)
fviz_pca_ind(res.pca,
col.ind = design_PCA,
repel = TRUE,
palette = "Set1",
ggrepel.max.overlaps = 1) + labs(title ="PCA normalisad count (DESeq2)",)
Heatmap d’expression de la région autour de POLR3G
ComBat_norm = norm_list_ComBat_melted$DESeq2
# VARIABLES
TNBC = c("MDAMB231","MDAMB436","HCC1937", "HCC1806",
"HCC1954","MDAMB468")
of_interest = c(
"POU5F1" ="ENSG00000204531",
"MYC" = "ENSG00000136997",
"POLR3GL"="ENSG00000121851",
"POLR3G" ="ENSG00000113356",
"MEF2C" = "ENSG00000081189",
"CETN3" = "ENSG00000153140",
"MBLAC2" ="ENSG00000176055",
"LYSMD3" ="ENSG00000176018",
"ADGRV1" ="ENSG00000164199")
of_interest = data_frame("ENS_id" = of_interest, "gene_name" = names(of_interest))
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Filtering matrix with only gene of interests
tmp = filter(ComBat_norm, grepl(paste(of_interest$ENS_id,collapse="|"),ComBat_norm$Var1))
tmp$copie = substr(tmp$Var1,1,15)
# merge dataframes to correctly match ENS_ID with gene name
to_heatmap = select(merge(tmp,of_interest, by.x = "copie", by.y = "ENS_id"),
gene_name,Var2 ,value)
colnames(to_heatmap) = c("gene_name","cell_line","expression")
# melt ==> classical matrix
heat = reshape2::dcast(to_heatmap, gene_name~cell_line)
## Using expression as value column: use value.var to override.
row.names(heat) = heat$gene_name
heat$gene_name = NULL
# Median centering
heat_norm = t(apply(heat, 1, MedianCentering))
# Create color scale
scale = get_color_scale_pheatmap(heat_norm,ncolors = 50,threshold = 0.75)
annotation_heatmap = data.frame("data_source" = design, row.names = cell_lines)
color = c("red","blue")
names(color) = c("CCLE","RoadMap")
annotation_heatmap$state = ifelse(
grepl(paste(TNBC,collapse="|") ,row.names(annotation_heatmap)),
yes = "TNBC", no = "other")
ComBat_heatmap = data.frame(t(heat_norm))
ComBat_heatmap$colors = ifelse(grepl(paste(TNBC,collapse="|") ,row.names(ComBat_heatmap)),
yes = "red", no = "darkgrey")
ComBat_heatmap = ComBat_heatmap[order(-ComBat_heatmap$POLR3G) , ]
ComBat_heatmap = ComBat_heatmap[,c("CETN3","MBLAC2","LYSMD3","MYC","POLR3GL","ADGRV1","MEF2C","POLR3G","colors")]
map = pheatmap(select(ComBat_heatmap,-colors),
color = scale$colorscale,
breaks = scale$breaks,
cluster_rows = F, cluster_cols = F,
angle_col = 0,
fontsize_row = 5,
main = "Expression heatmap : ComBat batch correction ")
#tmp = ComBat_heatmap
map$gtable$grobs[[4]]$gp=gpar(col= ComBat_heatmap$colors, fontsize = 5)
print(map)
limma_norm = norm_list_limma_melted$DESeq2
tmp = filter(limma_norm, grepl(paste(of_interest$ENS_id,collapse="|"),limma_norm$Var1))
tmp$copie = substr(tmp$Var1,1,15)
to_heatmap = select(merge(tmp,of_interest, by.x = "copie", by.y = "ENS_id"),
gene_name,Var2 ,value)
colnames(to_heatmap) = c("gene_name","cell_line","expression")
heat = reshape2::dcast(to_heatmap, gene_name~cell_line)
## Using expression as value column: use value.var to override.
row.names(heat) = heat$gene_name
heat$gene_name = NULL
heat_norm = t(apply(heat, 1, MedianCentering))
scale = get_color_scale_pheatmap(heat_norm,ncolors = 50,threshold = 0.75)
annotation_heatmap = data.frame("data_source" = design, row.names = cell_lines)
color = c("red","blue")
names(color) = c("CCLE","RoadMap")
annotation_heatmap$state = ifelse(
grepl(paste(TNBC,collapse="|") ,row.names(annotation_heatmap)),
yes = "TNBC", no = "other")
limma_heatmap = data.frame(t(heat_norm))
limma_heatmap$colors = ifelse(grepl(paste(TNBC,collapse="|") ,row.names(limma_heatmap)),
yes = "red", no = "darkgrey")
limma_heatmap = limma_heatmap[order(-limma_heatmap$POLR3G) , ]
limma_heatmap = limma_heatmap[,c("CETN3","MBLAC2","LYSMD3","MYC","POLR3GL","ADGRV1","MEF2C","POLR3G","colors")]
map = pheatmap(select(limma_heatmap,-colors),
color = scale$colorscale,
breaks = scale$breaks,
cluster_rows = F, cluster_cols = F,
angle_col = 0,
fontsize_row = 5,
main = "Expression heatmap : limma batch correction ")
#tmp = ComBat_heatmap
map$gtable$grobs[[4]]$gp=gpar(col= limma_heatmap$colors, fontsize = 5)
print(map)
